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A numerical technique, alternative to Grad's well-known ADM method has 
been proposed to deal with the slow adiabatic evolution of a toroidal plasma 
with flow. The equilibrium problem with prescribed adiabatic constraints may 
be solved by simultaneous calculations of flux surface geometry and original 
profile functions. Implications for the problem of bifurcation due to nonlinearity 
of the governing equations are discussed. In the case of field-aligned sub- 
Alfvénic flow the system is in the second elliptic regime if 8 < A?/(1— A?) at the 
magnetic axis, where A is the Mach Alfvén number of the flow. Super-Alfvénic 
flows do not satisfy the local firehose stability criterion. 


1. Introduction 


An axisymmetric ideal MHD equilibrium state with flow can be defined if one 
specifies five arbitrary functions of the poloidal magnetic flux y, which are 
denoted by F, H, w, ® and S (Zelazny et al. 1991). The function H is an analog 
of the Bernoulli function in fluid dynamics, S is the entropy of the plasma, and 
F determines the toroidal magnetic field, modified by flow. The functions ® 
and w are related to the field-aligned and toroidal flow velocities respectively. 
The latter is driven by the electrostatic field. These quantities are not observed 
directly in à physical experiment, although some of them have clear physical 
meaning (o is the frequency of rigid toroidal rotation of a magnetic surface). In 
this paper we present a formulation of the problem of toroidal plasma flow 
equilibrium in terms of five other functions, which turn out to be adiabatic 
constraints of slow plasma evolution. We consider a plasma in a toroidal device 
with nested magnetic surfaces, possibly surrounded by a vacuum region. The 
plasma edge is assumed to be a magnetic surface. In $2 we give a summary of 
the basic equations, followed in 83 by a presentation of the transformation 
procedure from the original magnetic flux functions to the physically relevant 
flux functions. We then, in $4, formulate appropriate boundary conditions, 
since the problem is now posed in flux co-ordinates with the magnetic axis 
as a singular origin of the system. The implications for the problem of 
bifurcation due to nonlinearity of the equations are briefly discussed in $5. A 
successful numerical simulation of a particular problem will be reported in a 
separate paper. Our approach is complementary to the Zehrfeld and Green 
concept of reference line (Zehrfeld & Green 1972; Semenzato, Gruber & Zehrfeld 
1984). 
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2. Governing equations 


Let us briefly summarize the equations that describe a stationary flow 
equilibrium. The ideal stationary MHD equations are (Kerner & Tokuda 1987) 


V.(pv) ^ 0, (1) 

p(v.V) v - Vp = (Vx B)x B, (2) 
Vx(vxB)- 0, (3) 

V.B — 0, (4) 

(v.V)S — 0, (b) 


where p is the mass density, v is the flow velocity, B is the magnetic induction, 
p is the plasma pressure and S is the entropy per unit mass. In the presence of 
an axial symmetry, with co-ordinates (R, Z, €) centred at the main axis of the 
torus, (4), (1) and (5) imply 


B = Vex Vy -IVE, 
v= P E Ral) VC, 


8 = S(y), 


where ® and w are surface functions representing the field-aligned and toroidal 
flows respectively, and y is the poloidal flux of the magnetic induction B. 
The V6 component of (2) gives 


_ F(y) - 99E? 


Hs 1—9?/p ' (6) 


where F is a surface function. 
Integration of the equation of motion (2) along B gives the Bernoulli 
equation: 
2p2 
HU) = rs Mu ey, (7 


where H is an arbitrary function of the poloidal magnetic flux only, and T is 
the plasma temperature. The Bernoulli equation serves to define the density p 
in terms of (R, y, |Vy/]). 

Finally, projecting (2) along Vy, we get the extended Grad-Shafranov- 
Schlüter (EGSS) equation (Zehrfeld & Green 1972; Hameiri 1983; Kerner & 
Tokuda 1987; Zelazny et al. 1991): 


p? Vy - IF’ , , Seam ^ 
v4( 2) --% —pH' t pTS' — v.B' —Lw', (8) 
where 
v.B = (pR?) (6 |Vy |* -- IL), (9) 
L = Io + R?po, (10) 


and a prime denotes differentiation with respect to y. The equation of state 


p =(y—1)pT = exp (5) p, (11) 
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has been used throughout, with the adiabatic constant y, usually taken to 
be 3. 

In this manner we have two equations, for y and p, in terms of the five 
surface functions: F, H, S, ® and w. 

Now, we shall consider the evolution of such an equilibrium due to an 
increase of the magnetic flux in the vacuum region. The time scale of the 
evolution is assumed to be much longer than the Alfvén wave transit time and 
much shorter than the diffusive transport time scale. The former means that the 
plasma evolves through equilibrium states, the latter that the evolution is 
subject to some adiabatic constraints, namely that the following quantities are 
conserved during the evolution of the system: 


cross-helicity o (y) — fv-Ba, 
toroidal momentum L(y) = fi dV, 

; : I 5 
toroidal magnetic field flux X (y) = Í 5;dV, 
mass My) = fear. 
thermal energy E(w) = fora, 


"where integration is over a volume enclosed by a given magnetic surface 
V = const. These constraints determine the right-hand side of the equilibrium 
equation (8). The derivatives of these functions with respect to y are 


val np. n 


dV dV dV 
(a, fag 
s = m (3p) d = (5 | 
where 
usw Dr db» k= (xs). 


m=<py, e=<pT), 


V(y) is the volume (divided by 27) of the region bounded by the surface 
V = const, and 


(y= [maag 


is a surface averaging operator, with 6 denoting the angular poloidal co- 
ordinate within the magnetic surface, and 


dV=RGdyd0, G —0,R0,Z—0,R0,Z. 
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In what follows we limit the problem to the case of the field-aligned flow, w = 0. 
A pure toroidal rotation is only significant if it is near-sonic (Semenzato et al. 
1985). 

Taking surface averages of (6), (7) and (9)-(11), we obtain the transformation 


formulae 
ye = Hm—}®w, from the Bernoulli equation (7), 


k=FF,, from (6), 
l=FOF,, from (10), 
e=exp(S)F,, from (11), 

w = O(F,+F?F,), from (9), 


^T cera) ^ cen) B= Ger) 


1 <p”? 
F, =( — F = ; 
^ aces?) aim 
The EGSS equation (if expressed in terms of the adiabatic invariants W (y), 
L(y), etc.) can now be written in the form 


where 


on a 5) 
A v= F(R. SS j (12) 
where A* is the extended Grad-Shafranov-Schlüter differential operator and 
the right-hand side is not given as an explicit function of y but as a functional 
4 (R,V,y,dy/dV,d*v/dV*) with V the volume (divided by 27) inside a 
magnetic surface y = const and the magnetic flux y being considered as a 
function of V. 
One can also surface-average the EGSS equation (8), obtaining 


1 d , * X "OE LN r ^ / 5] 
ral" (( 3 r2 |= kF’ —mH' + eS’ —wd' — lw’. (13) 


This is a dual form of the EGSS equation in the sense that the functions k, m, 
w, | and e are now given, while the functions F, H, ®, w and S are to be 
determined. 

Equations like (12) are called generalized or queer differential equations. 
Considering static equilibrium, Grad (see Grad, Hu & Stevens 1975) proposed 
a simple iterative technique called the alternating dimension method (ADM). 
With this method, one iterates between a two-dimensional problem of geometry 
calculation and a one-dimensional problem of profile calculation. For the 
geometry calculation the right-hand side of the Grad-Shafranov-Schlüter 
equation is considered given (from the previous iteration) and the magnetic 
surfaces y(r) are determined from the GSS elliptic equation. Then profiles are 
determined by solving the equations for profile transformation. This method 
has been applied to the problem of a flux-conserving tokamak (FCT) (Dory & 
Peng 1977). Hameiri (1983) generalized Grad's technique to the case of flow 
equilibrium. In this paper we present a method that is able to solve the 
problem without iteration between profile and geometry calculations. 
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3. Methodology 


It was proved by Baransky (1987) that the set of equations (1)-(5), in the 
presence of nested magnetic surfaces, has the property of self-similarity with 
respect to the entropy distribution, S(y). The other variables in the problem 
may be redefined so as to absorb this function, thereby eliminating it from the 
problem. A single solution of this transformed set of equations represents an 
infinite set of solutions to the original set of equations. This allows us to adopt, 
without loss of generality, an arbitrary entropy profile, for instance S = 0, i.e. 
exp(S) 21 and S'z0, which substantially simplifies the transformation 
formulae. We obtain 


k 
F TRES (14) 
P= His (15) 
EJ 
H = m` [yF, 39 (P,  F;F,)]. (16) 


Equation (15) is in fact an integral equation of the form 


of” RG dO -i Gdo 
o 1—-9*/p kJ, R(L—O*/p)’ 


with J and k given, and ® to be determined. 
Instead of solving (14)-(16), we differentiate them with respect to y to get a 
system of ODEs for the functions H, ® and F: 


k—FF, 
= 17 
F ROC (47) 
, F-9( SF FIF, 
® es (18) 
gi — [09 0, Py) H9", + 2FF F, + PF) + yF] Hm’ (i) 
—————O ieee 


As (17)-(19) for the original profiles are differential equations with respect to 
the poloidal magnetic flux y, it is quite natural, and to some extent necessary, 
to formulate the whole problem in a flux surface co-ordinate system (0, 0), where 
ô labels magnetic surfaces (i.e. V = y(ô)), ô = 0 at the magnetic axis and ô = 1 
at the plasma surface; 'is an angular co-ordinate within magnetic surfaces. 

In flux co-ordinates (6, 0) the EGSS and Bernoulli equations have the 
following forms: 


1 ô 2a ( P | did Al -z | 
koel p uL *gR 06 FAN p yi 


1 
&pl s (o Duy Je! = 0, (20) 


H= ju 24 72 o y y-1 91 
Val + Rp oe p > ( ) 
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where 

g—R,Z,—H,Z, ga =R +Z Js = — (Ry Ry t+ Z, Zo). 
The single subscripts denote differentiation with respect to the respective co- 
ordinates. 

Because of the intrinsic nonlinearity of the EGSS equation and strong 
coupling to the Bernoulli law (which is also nonlinear), approximate or 
numerical techniques need to be invoked to solve the governing equations 
except for some special examples of static equilibrium or unidirectional (in the 
toroidal direction) flow (Laing, Roberts & Whipple 1959; Infeld 1971; Maschke 
& Perrin 1980; Agim & Tataronis 1985). We adopt the moment method of Lao, 
Hirshmann & Wieland (1981) generalized to the case of flow equilibrium by 
Zelazny et al. (1991). The first step is to expand the co-ordinates (R, Z) in Fourier 
series in 0: 


R(8,0) = R, +acos6+..., (22) 
Z(0,0) = EasinO- ..., (23) 
where ‘...’ denote higher harmonics. 


In order to eliminate the derivative 0p/08 from (20), it is necessary to 
differentiate the Bernoulli equation (21) with respect to 6. Then, substituting 
(22) and (23) into (20), we obtain an ordinary differential equation of the form 

d'a ,d*E @R ,4,0E 5 
Basta TE —— + Ea cos 0 is 9 — a? cos rr RE =). (24) 
where Y comprises derivatives of lower order. The moments of (24) with respect 
to the basis (1, cos, cos? 6, ...) are 


2 ~ 
Hoe = 27,4... (25) 
2 ~ ~ 
Bate 49,7 ds... (26) 
2 ~ ~ 
ae a 41-81 (27) 
where 
Peo (fae Fo) life PSs Foro 
9 on ec Tee o o TR Bg hom ; 


The Bernoulli equation (21) serves to determine the density p in terms of 


(R, y», IVY). 

The above equations together with (17)-(19) constitute the set of coupled 
ordinary differential equations subject to the boundary conditions. The 
appropriate boundary conditions are derived in the next section. 


4. Boundary conditions 
We adopt the following labelling of the flux surfaces: 


y(ò) € V max(1 =O"); 


Toroidal equilibrium 135 


where Ymax is the value of the magnetic flux at the magnetic axis (ô = 0), 
whereas at the plasma edge (ô= 1) the boundary condition (y = 0) is 
automatically satisfied. We require the plasma boundary to be simultaneously 
an isobar and a magnetic surface, and therefore the flow must vanish there. The 
value Ymax of V at the magnetic axis has to be determined together with the 
solution. Thus the MHD equilibrium problem can be considered a generalization 
of an eigenvalue problem rather than a standard boundary-value problem. 
Because of this, the appropriate boundary-value problem has to be over- 
determined. For profile functions the boundary conditions are 


F(B = 0) -5-1u[ -2), 


= E = Y -1 1,2 
H(à = 0) = H, y-i^ T$. 


Here Jy, po, vo and By, are respectively the poloidal current, mass density, 
toroidal rotation velocity and toroidal magnetic field at the magnetic axis, 
where J, — B4,,E,, and Ro, is the location of the magnetic axis. These 
parameters together with the respective profile functions k, | and m constitute 
input to the problem. 

For the EGSS equation in the form of the ODE system (25)-(27) we 
formulate a two-point boundary value problem. 

(i) At the plasma boundary, à = 1, we fix the shape of the plasma surface: 


Ry = Ho, a = 0, E-E,.., 


where Ro» aj, Ep, ... are subject to slow evolution due to an increase in the 
magnetic flux in the vacuum region. 

(ii) At the magnetic axis, 0 — 0, we fix conditions that follow from 
requirement that our solution be analytic at the magnetic axis. We then assume 
the following expansions for functions Ry, a, E, ... and profile functions k, l, 
and m: 

Ry = Ro tRyt..., 
a =at... 


E =E, +E, +... 


bg ME o ee 
Roo 


l= htl ðt... lo = Po Vo Poo 
m = pg p3 0 +... 
The expansions of the profile functions k, l and m must be polynomials in even 
powers of ô to avoid singularities in the derivatives of the various physical 
quantities. This implies also that R, — 0 and E, — 0. Details are given in 
Zelazny et al. (1991). In deriving the last three formulae, we used the fact that, 


with the assumed expansions of the functions R,, a and E, the surface-averaging 
operator tends to the identity operator in the limit ó— 0. 
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From the surface-averaged EGSS equation (13), we have at the magnetic axis 


E -Se(i i) 
$-0 Roo p) Eo ` 


This is an extra boundary condition required for the determination of Ymax- 

Substituting (14)-(16), after some rearrangement, the following form of the 
term in parentheses on the left-hand side is found, with all derivatives taken at 
ó 0: 


d d? lo ko d 
a (t ap, gg o) (28) 


d? @ Ik, d 
koga” + Po gga ts 


d? 
- li — 


2 (9 9 2 d 
TES, elc ky Ro, | 1— 2) an- — 2k, 
l @?\ a? 
+2 |». -l (i -X a 
Do Boo po dò? * 


1 ly 2 s Q2 2 d? 
talaka) pr 7) aet 
OE M n 
2k, Ri då? * 
For not too high velocities, A? < 1, where A = [®| pt is the Mach Alfvén number 
of the flow, (28) takes the much simpler form 
sut AV ax E? +1 


d /k 
a? 
t [ass or "a so Re B 


With the aid of Mathematica (Wolfram 1991), the second derivative d*F,/d? can 
be calculated: 


+ 


dF,  a2,—8Ry Roe 
as (20) 


Standard numerical methods (e.g. the shooting method), can now be used to 
solve the boundary-value problem for the set of ODEs (17)-(19) and (25)-(27). 
The quantities Ymax E, and Roo may be chosen as shooting parameters, while 
a; is calculated from (28), which, in view of (29), is a biquadratic equation of 


the form 
kat, 20,03, — C, — 0, 


where 
C, = Al Roo Roz + ko ka Roo - Yio Pb’ 
o o dta LES 
TUR OEC 
so that 


E 4T (C1 C, Am 
E ERE - LOS 
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Following Lao et al. (1982), we take for R, and E, small negative values of the 
order of 107+. 

Our methodology differs from Grad’s alternating dimensions method (ADM) 
(see Grad et al. 1975) in the sense that there is no iteration between geometry 
calculation and profile calculation, but instead profiles are calculated 
simultaneously with geometry. 


5. Discussion 


It is well known (Zelazny et al. 1991) that, owing to the nonlinearity of the 
Bernoulli equation, the EGSS problem stated in terms of the original profiles H, 
®, w, S and F may have up to four distinct solutions, which are distinguished 
as first elliptic or first hyperbolic, second elliptic, third elliptic and second 
hyperbolic. Namely, if A? is the square of the Mach Alfvén number of the field- 
aligned flow and F = £/(#+1) is the speed of the trailing cusp of the slow 
wave, with f = yp/B?, then for 0 € A? <T the system is in the first elliptic 
regime. For T < A? < A?, where A, is the Mach Alfvén number of the slow 
magneto-acoustic wave with respect to the poloidal Alfvén speed, the system is 
in the first hyperbolic regime. The second and third elliptic regimes correspond 
to A} < A? « 1 and 1 « A* « Aj respectively, where A, is the Mach Alfvén 
number of the fast magneto—acoustic wave. Finally, for A? > A}, we obtain the 
second hyperbolic regime (see figure 1). The solutions for p from the Bernoulli 
law correspond to the intersections of the curve shown in figure 1 with the 
horizontal line H — const. The solution chosen accordingly affects the nature of 
the differential operator in the EGSS equation. 

We should now like to interpret the above classification of the EGSS 
equation, shown in figure 1, in terms of the physically meaningful parameters 
A? and f. With this purpose in mind, let us recall that the transition from the 
first hyperbolic regime to the second elliptic regime takes place at values of 
the parameters satisfying: 


D? 1-5—[(L4- B? —48B2/ B") 


p 2B? /B? 
At the magnetic axis (B, = 0) this reduces to 
9$ f, 
Po Pott 


where f, = yp,/ B,. In terms of the parameters of $4, this condition is 


2 1 
40 | 42 = 1, A, = xA 
0 TO 
or, more generally (for B, > 0, i.e. also in the off-axis region), 
Bp As) 4e 
(i- 24 \G+4 =1, 
where A is the Mach Alfvén number of the flow, 
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ra 4 


Figure 1. Plot of the function given by the right-hand side of the Bernoulli law versus the 
square of the Mach Alfvén number, A*: eI, eII and eITI indicate the first, second and third 
elliptic regimes; hI and hII indicate the first and second hyperbolic regimes. 


So, if 8, < A2/(1 — A3), we have the second elliptic solution at the magnetic axis, 
and the first elliptic solution otherwise. Additionally, in order to avoid entering 
the first hyperbolic regime (for which the boundary-value problem should be 
formulated in a quite different manner), the following conditions should be 
fulfilled in the off-axis region: 


2 2 
p« ( Ba)“, for the second elliptic regime, (30) 
A? ae : 
p? 14 for the first elliptic regime. (31) 


Second hyperbolic and third elliptic regimes are reached for super-Alfvénic 
plasma flow velocities, A? > 1. These solutions do not meet the local firehose 
stability criterion o = 1+ (p, —p,)/B’ > 0 if p, is interpreted as p pv? and p, 
as p. The condition ø > 0 corresponds to the requirement v < v,, which means 
sub-Alfvénic speeds for the field-aligned flow. We conclude that if flows are 
present, the second elliptic solution occurs for a low-£ plasma (cf. (30)), and the 
first elliptic solution for a high-f plasma (cf. (31)). Note that second elliptic 
solutions are additionally distinguished by the local criterion for mirror-mode 
stability, 7 = 1--2dp, /dB* >0 (where the derivative should be calculated 
along B; see e.g. Fowler 1981), which, in terms of f and A?, takes the form 
p<A?=>7>1, 
2 
A*«f« a =>7<0. 


For the first elliptic branch, (31) implies 0 < 7 < 1. This means that the first 
elliptic solution and the second elliptic low-f solution (f < 4?) are stable 
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0 1 2 3 
Square of Alfvén Mach number, A? 


Figure 2. Regimes of ellipticity and hyperbolicity of the EGSS equation: A, first elliptic; B, 
first hyperbolic; C, second elliptic; D, second elliptic; E, third elliptic; F, second hyperbolic. 


against local mirror modes (see regions A and D in figure 2). This discussion is 
summarized in figure 2, which is drawn for B2/B* = j. 

These results suggest that, if flows are present, it may be easier to find a 
plasma in the second elliptic regime than in the first elliptic regime, even though 
it is the former that may be reached continuously from the static equilibrium. 
It was shown by Zelazny et al. (1992, 1993) that solutions of the second elliptic 
branch have density axis markedly shifted inwards with respect to the 
magnetic axis and that they correspond to a more elongated plasma cross- 
section. This has to be taken into account in the equilibrium reconstruction 
from magnetic measurements, which is done, to the best of our knowledge, in 
all cases on the basis of static equilibrium with density and magnetic axes 
coinciding. 


The authors wish to thank Drs R. Stankiewicz and S. Potempski for 
beneficial discussions. This work has been supported by Grant 201939101 from 
KBN in Poland. 
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